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We present an algorithm to compute the number of so- 
lutions of the (constrained) number partitioning problem. A 
concrete implementation of the algorithm on an Ising-type 
quantum computer is given. 

PACS numbers: 03.67.Lx,89.70.+c 



I. INTRODUCTION 

The discovery of quantum algorithms (QA's) that, 
when executed on a quantum computer (QC), give signif- 
icant speedup over their classical counterparts has 
given strong impetus to recent developments in the field 
of quantum computation. In theory an ideal QC is a 
universal computer. This means that for a given prob- 
lem there exists an algorithm to solve this problem on 
a QC. The fact that a QC is a universal computer does 
not tell us the algorithm itself, nor the computational 
resources that are required to solve the problem. 

In general it is not easy to construct algorithms for 
an ideal QC. In particular, algorithms that involve many 
conditional jumps (e.g. IF-THEN-ELSE statements) are 
difficult to implement. In essence this is because test- 
ing for a condition requires a measurement during the 
execution of the program. It is therefore of interest to 
see how a QC can solve problems that a conventional 
computer solves by performing many conditional jumps. 
The purpose of this paper is to present a new QA for one 
such problem of combinatorial optimization: The (con- 
strained) number partitioning problem. 



II. NUMBER PARTITIONING 

The number partitioning problem (NPP) is defined as 
follows l^-U]: Does there exist a partitioning of the set 
A = {ai, . . . , a n } of n positive integers cij into two dis- 



same footing by asking if there exists a partition such 
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< A where A = 1 (0) if B 



joint sets A\ and Ai = A — A\ such that £ a 



EA 2 



? The answer to this question is trivially no 
Sa,eA a 3> is odd - More g en " 



if the sum of all a.,-, B 



erally, the case of even or odd B can be treated on the 



is odd (even). 

For certain applications there may be additional con- 
straints on the partitioning of the set A. A common one 
is to fix the difference C between the number of elements 
in Ai and A 2 : C = £^ ei4l 1 - E Qj eA 2 1 - For instance, 
if C — we ask if there is a partitioning such that the 
number of elements in A\ and Ai are the same. 

For a given instance of A = {ai, . . . , a„}, we may en- 
code the whole problem using only n log 2 B bits. The 
NPP can be solved by dynamic programming, in a time 
bounded by a low order polynomial in nB M. As nB is 
not bounded by any polynomial of the input size n log 2 B, 
the dynamic programming algorithm does not solve the 
NPP with polynomial computational resources Q . 

Number partitioning is one of Carey and Johnson's six 
basic NP-complete problems Q. In practice, a problem 
is NP-complete if its solution requires computational re- 
sources that increase faster that any polynomial of the 
input size. Number partitioning is a key problem in the 
theory of computational complexity and has a number 
of important practical applications such as job schedul- 
ing, task distribution on multiprocessor machines, VLSI 
circuit design to name a few. 

The computation time to solve a NPP depends on the 
number of bits b = log 2 B needed to represent the in- 
tegers dj and B. Numerical simulations using random 
instances of A show that the solution time grows expo- 
nentially with n for n -C b and polynomially for n b 
|^-^|. For random instances A, the NPP can be mapped 
onto a hard problem of statistical mechanics, namely that 
of finding the ground state of an infinite-range Ising spin 
glass [pi)|-[l^|. The transition from the computationally 
"hard" (exponential) to "easy" (polynomial) has been re- 
lated to the phase transition in the statistical mechanical 
system jnjjij]. 

Although the transition between easy and hard prob- 
lems is important from conceptual point of view, it is 
good to keep in mind that most real-life problems are 
of the easy type jf|. For instance, if the a/s represent 
the weight of boxes that are to be distributed over several 
trucks, it is highly unlikely that the weight of these boxes 
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will vary between say 1kg and 2 32 kg, or that it is impor- 
tant to know the weights of the boxes with a precision of 
e.g. ten digits. 



III. QUANTUM ALGORITHM 
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1 — exp 


2ti"*(E"=i Oj^j - A) 




1 — exp 





(3) 



The potential power of a QC stems from the fact that a 
QC operates on superpositions of states 13|-ig[|. The in- 
terference of these states allows exponentially many com- 
putations to be done in parallel fll3|-p"9|| . A QA consists 
of a sequence of unitary transformations that change the 
state of the QC @^9[. Therefore to solve a NPP on 
a QC, we first have to develop an algorithm that does 
not contain conditional branches and can be expressed 
entirely in terms of unitary operations. 

A generic n-qubit QC can be modeled by a collection 
of n two-state systems, represented by n Pauli-spin ma- 
trices {(?!,...,(?„} [l^-|ist|. The two eigenstates of cr? will 
be denoted by | f) ■ and | J.) •, corresponding to the states 
\0)j and of the j-th qubit respectively. The eigen- 
values corresponding to | |)j and | J,) ■ are Sj = +1 and 
Sj = — 1. They can be used to represent a partitioning of 



A: We assign a,j to A\ (A 2 ) if Sj 
can find a set {Si, . 



If we 



jSj — 



-1 (Sj = -1) 
, S n } such that A — Ej=i a 
we have found one solution of the NPP. 

It is known that the most simple class of spin systems, 
i.e. those involving interactions of the Ising type only, 
can be used to build universal QC's |l4|Jl^ , |20| . For our 
purposes it is, at this stage, sufficient to consider a sys- 
tem of n non-interacting Ising spins. The Hamiltonian is 
given by 



As | A — Ej=i a jSj | < M for any choice of {Sj}, the sum 
over to in (2) will be zero unless A — Y^j=i a j^j = 0> 
in which case the configuration {Si, . . . , S n } is a solution 
of the NPP (note that there can be exponentially many 
solutions, for instance in the case that all the a/s are 
equal). Performing the sum over all spin configurations 
as indicated in (^), it follows immediately that n s is the 
number of solutions of the NPP. Note that (2) gives the 
number of solutions of a NPP, which is more than just a 
yes or no answer to the question if a partition of A exists 

Formally expression (H) is the density of states at zero 
energy of the physical system described by Hamiltonian 
(jl]). Elsewhere we have shown that for a large class of 
models H, the density of states can be calculated effi- 
ciently on a QC p2fl . The algorithm presented below, 
although related to the one described in p2|, is specifi- 
cally tuned to solve NPP's. 

The equivalence of (2) and the solution of the NPP 
can also be demonstrated by explicit calculation of the 
trace over all spin configurations. This is easy because 
the spins do not interact. The result is 
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where the a^'s represent external fields acting on the 
spins. From it follows directly that an eigenstate of 
H with energy zero corresponds to a solution of the NPP. 
We will use Hamiltonian (Q) to define the time evolution 
of the QC, i.e. the QA that solves NPP's. 

The second key to the construction of the quantum 
algorithm is the observation that the number of solutions 
n s of a NPP is given by 
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M-l 



m— 



(2) 



where M = B + A + 1 and Tr U denotes the trace of the 
matrix U . Using the representation that diagonalizes 
the spin operators cr? we find 



For A = and in the limit M 
where 



oo we have n, = 2 n I s 



1 

2^ 



cos(ai6') . . . cos(a r 



(5) 



The question whether I s = or not is known to be equiv- 
alent to the (non-)existence of a solution of a NPP p.p3[. 

If n s > we can find a partitioning in the following 
manner. Assume we already know the values of the first 
< I — 1 < n spins. We make a guess for Si and compute 
ni° = M _1 J2m=0 * r e- 2mmH/M where the use of the 
symbol tr instead of Tr indicates that in calculating the 
trace, the values of the variables Si, . . . , Si are fixed. If 
ni > our guess for Si was correct, if not we reverse 
Si- Then we increase I by one and repeat the procedure. 

The algorithm outlined above is easily generalized to 
handle constraints. Introducing another Hamiltonian 
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the number of solutions n s (C) of the constrained NPP is 
given by 
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2-nimH / M -2irikH' / K 
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where K = n + \C\ + 1. Repeating the same steps as 
above we find that the sum over k yields zero unless 
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corresponding to m) reads 
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The procedure to find a partitioning itself is the same as 
in the unconstrainted case. 

The algorithms defined by (0),(|) and (§),(§) solve 
NPP's and constrained NPP's without recourse to dy- 
namic programming. This follows directly from explicit 
expressions (|J) and (||) . On a conventional computer the 
computation time required is bounded by nM (or nMK 
for the constrained case). Hence also these algorithms 
do not solve the (constrained) NPP in polynomial time 
(space). The conceptual difference between these algo- 
rithms and the dynamic-programming approach is that 
the former directly compute the number of solutions of 
the NPP whereas the latter performs a search for a solu- 
tion of the NPP. 

As we now show, (^) (or ([?])) is a convenient starting 
point to construct a QA that solves the (constrained) 
NPP on a QC. As will be clear from the discussion below, 
it will be sufficient to concentrate on (Q), the algorithm 
for (0) is almost identical. 

The first step is to introduce a "number operator" X 
with eigenstates \x), X\x) — x\x), x — 0, 1,...,M — 
1. We modify the Hamiltonian that governs the time- 
evolution of the QC as follows: 



H = AX a 3 a]X. 

3=1 



(9) 



Calculating the trace in the basis that diagonalizes 
erf ... cr^ and X we find that n s = M^Tr e ~ 27Tiil / M . 
Because H is diagonal in this basis Tr e - 2 * lH / M i s pro- 
portional to one matrix element, namely 

e -2iriH/M _ 

2 n M(U 1 . . . UnUxle-^/Mp! . . . U n U x ), (10) 



where \Uj) = (| t)j + I i)j)/V2 is the uniform super- 
position of the spin up and down state of spin j, and 
\U X ) = (|0> + |1) + . . .+ \M - 1))/VM is the uniform su- 
perposition of all the eigenstates of the number operator 
X. To derive expression (|l0|) we made use of 



^\U j )=cos(a)\U j )-isin(a)\U j ), 



(11) 



and {Ui\Ui 



From 



> = 0, where \Uj 
it follows that 



(I T), 



i) 3 )/V2. 



2 n {U t . . . UMe-^/^Ui . . . U n U x ). (12) 



As a QC can compute e~ ttH \t(j) with one operation (for 
arbitrary input \ip}) p3[ , ( p^ ) shows that once the QC is 
in the state of uniform superposition \Ui . . . UnU x ), one 
time-evolution step of the QC will solve the NPP. 

The initial state | |, . . . , |, x = 0) can be transformed 
into the state of uniform superposition \Ui . . . UnU x ) by 
the standard procedure p6 17|: The states \U) ■ can be 



obtained from the initial state 
spin j about the y-axis, i.e. 



| t)j by a rotation of the 

\U) j = e- l ™V A \ T) j for 
j = 1, . . . , n. On an Ising-type QC the states can be 
implemented by adding new two-state systems. We de- 
note the corresponding Pauli-spin operators and eigen- 
values by jlp and s p respectively. We use these spins 
to represent x = JZt—i 2 l ~ 2 {\ — s{) in binary form. As 
< x < M the number of spins p required to represent 
x is the smallest integer p for which M < 2 P . Using this 
binary representation for \x), the uniform superposition 
\U X ) can be obtained by p rotations of the initial state: 



\U X ) = e —^/4| T)l 



-«TM?/4| 



(13) 



where ® denotes the direct product operation. The sys- 
tem now comprises n+p spins and its Hamiltonian reads 

p n n p 

« = -EE J i,i°*iti - E h°i - E + d, (w) 



i=i j=i 



i=i 



bj = aj(2 p -l)/2, q = A2 



1-2 



and 



where Jjj — 
d = A(2P-l)/2. 

The complete QA for computing n s , i.e. for solving 
NPP's, can be summarized as follows: The initial state 
of the QC (all spins up by convention) is transformed into 
the state of uniform superposition. This takes n + p one- 
qubit operations. Next the QC makes one time-evolution 
step exp(— inl-l/2 p ~ 1 ). The matrix element in ( fL2| ) is 
obtained by applying the inverse of the n+p rotations, 
followed by a projection onto the initial state. Clearly 
the total number of QC operations is only 2n + 2p + 1 
while the amount of memory used is C(log 2 M + log 2 n). 

The constrained NPP can be solved in the same way: 
Add qubits to represent the variable k in (Q) and repeat 
the steps that lead to (|lj). Note that once the uniform 
superposition has been prepared, the QA also solves the 
constrained NPP with one time-evolution step. 



IV. IMPLEMENTATION ON A QUANTUM 
COMPUTER EMULATOR 

For the purpose of demonstration we have imple- 
mented the QA that solves the unconstrained NPP on 
a Quantum Computer Emulator (QCE), a software tool 
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for simulating physical models of QC's p4|. A subtle 
point thereby is that ( |l2| ) is not directly measurable be- 
cause e ~ 27T iH/M j g no ^. a physical observable. However it 
is not difficult to express n s in terms of an expectation 
value of a physical observable. 

Let us write the number of solutions (7) as n s — 
2 n (0|$) where |$) = U^e-^W^U^) and U = 



T»/4 



ur aim is to 



replace the projection onto the initial state |0), a short- 
hand notation for the state with all spins up, by the mea- 
surement of some observable. This can be accomplished 
by introducing another spin k, initially in the state of 
spin up, to the system and flip this spin if the other n+p 
are all up, i.e. by performing an AND operation on the 
n+p spins. With V denoting the unitary transforma- 
tion that performs this AND operation, we have in the 
language of qubits instead of spins 



|*) = vU- 1 e-^ n/2P V|0) ® |0)„ 
= y[2-»n s |0)®|0) K + (...)®|0) f 
= 2-"n s |0)®|l) K + (...)®|0) K , 



(15) 



where |^) is an element of the direct product of the 
Hilbert spaces spanned by the n+p spins and the aux- 
iliary spin k. We use the abbreviation (...) to represent 
the sum of all states of the n+p spins that have at least 
one spin down. From (15) it immediately follows that 



2 n (*|(l- K z )/2|*> 1/2 . 



(16) 



It is well-known how to implement the AND operation 
on a QC pq ]. In our practical implementation |2(|, we 
have choosen to use a three-bit network, the Toffoli-gate, 
as a building block for realizing the AND operation on the 
n+p qubits pf| . By adding extra work qubits the com- 
plete network requires of the order of \og 2 (n+p) steps and 
n+p extra qubits to perform the AND operation. Clearly 
this does not change the polynomial time and space char- 
acter of the QA that solves NPP's. A block diagram of 
the complete quantum program is shown in Fig.l. We 
have implemented the QA on a 15-qubit QC and used it 
to solve the NPP's A = {1,2,3,4}, A = {1,1,1,4} and 
A = {2,2,2,4} (these examples are included in the soft- 
ware distribution |p6| ) . In the final state the expectation 
values of the 15-th qubit are 0.015625, 0.00390625, and 
respectively. The corresponding number of solutions is 
n s = 2, n s = 1 and n s = 0. Clearly the demonstration 
program correctly solves NPP problems. 



V. ALTERNATIVE IMPLEMENTATION 

The implementation described above has the same log- 
ical structure as other QA's p|,^,^9[: Prepare the QC in 
a state of uniform superposition, perform some unitary 
transformation to encode information and then apply a 
filter to extract the answer. We now show that there is 



another QA that solves the NPP but does not fit into 
this general scheme in that the first step is missing. 

Consider the time-dependent n-spin correlation func- 
tion 



C(t) = ($| e 



iH x t—z 



.a„e 



-iH x t 



(17) 



where H x — — Yjj=i a 3 a j The state |$) can be any n- 
spin state that is an eigenstate of of ... cr^, e.g. the state 
with all spins up. As the cj's are unitary operators, it is 
a simple matter to write down a QA that computes C(t) 
on a QC. Obviously C(t) is a physically observable quan- 
tity but it may require a rather complicated experimental 
setup to measure this n-spin correlation function. 

Substituting into ( |l7| ) the equation of motion for each 
spin, i.e. e 4ffx *erje~'~ 
find 



<r| cos(oji) — Oj sbx{ajt), we 



C(t) = <$| n;=i[l cos( aj t) - ia? sin( aj t)m 

n 

= Y\_cos(a k t). 

3 = 1 



(18) 



The Fourier transform of C(t) at zero frequency is di- 
rectly proportional to I s and hence to n s : 



t C(t)dt = 2- n n s 5(uj) + R(uj), (19) 



where the regular part R(to) is zero at cj = 0. From ( |l9| ) 
it is clear that the NPP has a solution if S(ui) shows a 
peak at zero frequency. Detection of the central peak 
in the dynamic correlation function S(tu) may require a 
very long observation time T. To distinghuish between 
n s = and n s — 1 the observation time T must be larger 
than 2™7r. 



VI. DISCUSSION 

The essence of the algorithm proposed above is that it 
expresses n s in terms of the density of states of a physi- 
cal system (Ising spins in our example) . Clearly this QA 
certainly has its weakness if n s is close to zero and n 
is large. Of the order of 2™ measurements of k z are re- 
quired to distinguish between n s = 0, and n s = 1. This 
is tantamount to random sampling. By formulating, as 
we did, the outcome of the calculation in terms of (an av- 
erage of) physical observables (see ( p^ ) or (|l^) instead of 
a (collapsed) state, this problem of efficiency is difficult 
to overlook [ ^7| . 

NP-completeness of the NPP refers to non-polynomial 
behavior as a function of nB (B being the input size). 
Our QA is polynomial in this respect. However, in the 
hard case (b 3> n) and for large n, to obtain a yes-or-no 
answer tremendous precision is required. In the absence 
of any information of the a^'s other than that they are 
positive integers, the range of n s extends from zero to 
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( n ™ 2 ) • Any algorithm that computes n s (on a conven- 
tional or QC) should be able to cover this range (oth- 
erwise it can never return the correct n s ). This implies 
that the whole computation has to be done with at least 
the same (high) precision. 

As the NPP example shows, a realistic assessment of 
the potential power of a QA should include a quantita- 
tive estimate of the precision and other computational 
resources (e.g. energy) that are required to obtain the 
correct answer. For our QA an estimate of the required 
precision follows from ([l6]) . Note that the alternative im- 
plementation yields a similar, physically equivalent, esti- 
mate for the observation time T. 

The range of numbers a physically realizable QC will 
be able to handle is directly related to the energy range 
in which the QC operates {—B to +B in the NPP case). 
Although not a problem of principle, the physics of QC 
hardware will definitely impose some constraints on the 
range of numbers. 

There are two other, potentially large, numbers in- 
volved in an NPP problem: The number of states N s = 
2™ and the number of computational units N cu (of mi- 
croscopic size) in the physical realization of the QC. 
There are two cases to consider. 1) Theoretical (com- 
puter science): We have to examine the worst case. Then 
N s ^> N cu so that our NPP algorithm has little merit. 
2) In practice: In numerical experiments n < 32 

(N s < 2 32 ) whereas for instance in NMR QC experiments 
N cu w 10 18 ^> N s psfi . Using a sufficiently large number 
of computational units and efficient detectors it should 
be possible to distinguish between n s = and n s = 1. 
Assuming that clever engineering can produce spin sys- 
tems such as dl4| ) , our Q A might be used to demonstrate 
that a physical QC can solve a non-trivial problem. 

This work is partially supported by the Dutch 'Sticht- 
ing Nationale Computer Faciliteiten' (NCF). 
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FIG. 1. Block diagram of the quantum algorithm that 
solves the number partitioning problem in polynomial time 
and space. In this example the first n — 4 qubits are used to 
represent the integers to be partitioned. The p — 4 qubits 5 to 
8 are used to determine the number of solutions of the num- 
ber partitioning problem. The remaining 7 qubits are used to 
relate n 3 to a physically measurable quantity: The expecta- 
tion value of the 15-th qubit. The unitary transformation U 
prepares the uniform superposition of the first 8 qubits, U is 
the inverse of U, and the combination of INVERT and AND 
gates sets the 15-th qubit to one if and only if the first eight 
qubits are all one. 
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